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We consider the fate of the Dirac points in the spectrum of a honeycomb optical lattice in the 
presence of a harmonic confining potential. By numerically solving the tight binding model we cal- 
culate the density of states, and find that the energy dependence can be understood from analytical 
arguments. In addition, we show that the density of states of the harmonically trapped lattice 
system can be understood by application of a local density approximation based on the density 
of states of the homogeneous lattice. The Dirac points are found to survive locally in the trap as 
evidenced by the local density of states. They furthermore give rise to a distinct spatial profile of a 
noninteracting Fermi gas. 



PACS numbers: 37.10. Jk,05.30.Fk 

I. INTRODUCTION 

Graphene is a carbon monolayer with a honeycomb 
crystal structure, which was only recently produced [1 . 
The band structure of graphene is intriguing in that the 
dispersion is linear in the vicinity of the Fermi energy. 
This makes the material a zero-gap semiconductor, with 
quasiparticles behaving as massless Dirac fermions, thus 
opening the possibility of studying quantum electrody- 
namics with electrons in a solid state system [2 . The 
existence of carriers described by the Dirac equation has 
been confirmed experimentally along with the demon- 
stration of an anomalous quantum Hall effect [S', W. 
The striking electronic properties of graphene makes it 
an interesting system not only for studying fundamen- 
tal physics, but also as a platform for device fabrica- 
tion [HE]. 

Building on the potential of graphene as a test bed 
for relativistic quantum theory, several theoretical pa- 
pers have pointed out that ultracold atoms in a hon- 
eycomb optical lattice could prove an attractive, alter- 
native system for simulating relativistic physics [THIS]. 
An optical lattice is a periodic potential, formed by in- 
terfering laser beams, in which atoms exhibit the same 
Bloch band physics as solid state electrons. But contrary 
to a solid state crystal both the depth and the geome- 
try of an optical lattice potential can be controlled by 
adjusting the intensity and configuration of the lasers. 
Hence an optical lattice provides a pristine environment 
for implementing condensed matter models, and probing 
many-body dynamics such as the superfluid to Mott in- 
sulator transition [TT, TB . In addition, from a quantum 
simulator point of view ultracold atoms posses the ad- 
vantageous qualities of controllable interactions (using a 
magnetic field tunable Feshbach resonance [16 ) and the 
possibility of mapping the rich internal state space of the 
atoms onto multiple spin degrees of freedom. By apply- 
ing additional light fields an artificial gauge field can be 



engineered [TT', TS^ . With non- Abelian gauge fields dif- 
ferent topological phases can be engineered [13 . Other 
schemes for producing a relativistic dispersion with opti- 
cal fields have also been proposed p!9H22 | 

However, in experiments the discrete translational 
symmetry of the optical lattice is broken by a trapping 
potential, which confines the atoms. This complicates the 
comparison with solid state phenomena, but it is often 
surmised that the confining potential is slowly varying, 
and that a local density approximation can therefore be 
used. 




FIG. 1: (Color online) Schematic view of the honeycomb lat- 
tice with a superposed harmonic trapping potential. The har- 
monic potential is centered on one of the lattice sites. 

In this paper we test the validity of this presumption 
and determine how the physics of graphene is modified 
in an inhomogeneous honeycomb optical lattice. Similar 
calculations have been done for cubic lattices in one and 
two dimensions [23H27] . In the present work we first give 
a brief review of how a honeycomb lattice potential can 
be generated in an experiment (for a longer discussion see 
e.g. [I2])- Then we consider the situation illustrated in 
Fig. [l] where the translational symmetry of the lattice is 
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broken by a parabolic offset of the site energies. We solve 
a tight binding model numerically for a finite system and 
characterize the spectrum by the density of states. The 
spectral features can be understood by a combination of 
analytical arguments and a local density approximation. 
We address the existence of Dirac particles in the inho- 
mogeneous lattice by plotting the local density of states 
and by calculating the density distribution of a noninter- 
acting Fermi gas in the combined lattice and harmonic 
potential. 

II. CONSTRUCTING A HONEYCOMB 
OPTICAL LATTICE 

A honeycomb optical lattice can be constructed by su- 
perposing three laser beams with wave vectors {i = 
1,2,3) of identical magnitude /c^ = '^tt/Xl lying in the 
x-y plane at 27r/3 angles with each other. If the three 
lasers have the same intensity and are linearly polarized 
in the z-direction, this gives rise to a lattice potential of 
the form [H [28] 

VL{r) = Vo [cos(bi • r) + cos(b2 • r) + cos((bi + bs) • r)] , 

(1) 

where bi = ks — ki,b2 = ki — k2 are the reciprocal 
lattice vectors. The lattice depth Vq depends on the in- 
tensity and the detuning of the lattice lasers, and here we 
only consider Vq > corresponding to a positive detun- 
ing, which produces the honeycomb structure of graphene 
with a spacing between nearest neighbor lattice sites of 
ao = 2Xl/ a/27. A negative laser detuning generates a tri- 
angular lattice potential. Using additional confinement 
along the z-axis an effectively two-dimensional system 
can be realized. 

The honeycomb optical lattice described above can be 
generalized in two straightforward ways. First, if the 
direction of polarization of the lasers is changed from 
perpendicular to the lattice plane to coplanar with the 
wave vectors of the beams, the resulting periodic light 
field is circularly polarized at the positions of the lat- 
tice minima, with the lattice sites forming an alternating 
hexagonal pattern of cr+ and cr~ polarizations [29 . In 
such a lattice atoms in different internal spin states will 
experience different light shifts [30 , and for atoms with 
spin projection |mi?| > the lattice potential becomes 
a periodic array of offset double wells ^1]. Secondly, if 
the laser intensities differ, an anisotropic honeycomb lat- 
tice is generated where the tunneling rates depend on 
direction. If the intensity imbalance is sufficiently large 
this induces a band gap in the single-particle spectrum, 
equivalent to the Dirac fermions acquiring mass [7[ [T^ . 

In the following we restrict our attention to the spin- 
independent, isotropic honeycomb lattice. However, the 
form of the lattice potential above assumes that the three 
lasers beams are plane waves, while in reality their cross 
sections have a gaussian intensity profile. This gives rise 
to an energy offset between different lattice wells, which 
for a sample much smaller than the beam widths can be 



approximated by a harmonic oscillator potential. In ex- 
periments an additional confining potential is often added 
intentionally to restrict the size of the cloud. Our moti- 
vation for this work is to investigate how the presence of 
such a spatially dependent energy offset between lattice 
sites affects the single-particle physics of the honeycomb 
lattice. 



III. TIGHT BINDING MODEL 

We consider a tight binding model with nearest neigh- 
bor tunneling and expand the Hamiltonian in terms of 
the localized (orthogonal) Wannier states of the first 
Bloch band, 

where \wj) is the Wannier state localized at lattice site 
j, and Tj is the distance of site j to the center of the 
trap, which has spring constant Hi. The sum in the first 
term is over nearest neighbor sites. The nearest neighbor 
tunneling amplitude between sites j and j' is defined as 

J=-{wj,\f + VL\wi) (3) 

where T is the kinetic energy operator in x?/-plane. Tun- 
neling to next-nearest neighbor sites is strongly sup- 
pressed. The tight binding model is illustrated in Fig. [l] 
For simplicity we take the center of the trap to coincide 
with one of the lattice sites. This restriction is easy to 
relax. 

We assume that the harmonic potential does not mod- 
ify the nearest neighbor tunneling rate. This approxi- 
mation is valid provided two neighboring wells (at dis- 
tances r and r -\- Sr from the trap center) are separated 
by a barrier Vb, which is much larger than the energy 
difference between their minima 3E(r). For r ^ ao the 
energy difference between neighboring points is 5E{r) < 
\Hi[{r + ao)^ — r^] ~ K.aQr. This defines an energy cutoff 
in our model, since 5E(r) increases with r. Thus high 
energy states with a wave function, which remains finite 
beyond a critical distance Vc = Vo//^ao, will not be rep- 
resented accurately in our model. Hence we are limited 
to consider energies E <^ Ec = \i^t1 — 3J, where — 3J 
is the lowest energy in the spectrum for a homogeneous 
lattice (see below). The relevant energy scale is set by 
the tunneling, and we thus require Ko^j J <C (Uq/J)^. If 
we introduce the characteristic length scale of the har- 
monic oscillator aosc = \fjj^ this criterion translates 
into ao/aosc <C Vb/J. 

The oscillator length scale is typically of the or- 
der of micrometers, while the lattice lasers have wave 
lengths of several hundred nm. Since Vq/J ~ 
(^o/^i?)'/^exp[1.582yV^], where Er = hy2mXl 
is the recoil energy of the lattice lasers [12], the condi- 
tion for the validity of the tight binding model is almost 
always satisfied for lattices deeper than about bEn. 
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In the numerical diagonalization we impose hard wah 
boundary conditions at Tmax = 60ao. This artificial re- 
striction leads to finite size effects, such as edge states, 
that may be interesting in their own right [TF. Below 
we also give analytic results, which apply for an infinite 
lattice. 



A. Homogeneous lattice dispersion 

We first give a brief review of the homogeneous lattice 
case where Hi = 0. For an infinite lattice the eigenstates 
of the tight binding Hamiltonian are Bloch waves with 
energies 



±J 



3 + 2 cos{\/3qyao) 



+4 cos ( -q^ao 1 cos I — %ao 



1/2 



(4) 



as a function of the quasimomentum q. The spectrum 
consists of a lower and an uppper band as depicted in 
Fig. |2] with a hexagonal first Brillouin zone. Near the 
six corners of the first Brillouin zone the two bands form 
opposing cones, which exactly touch at the corner points. 
Since each of the corners is shared equally between three 
adjoining Brillouin zones the first Brillouin zone contains 
two independent corner points at quasimomenta K and 
K'. Around these points the dispersion is linear: ^ 
±fi'UF|k| for q = K+k with |k| <C |K| and similarly in the 
vicinity of K^ Since this corresponds to the dispersion 
of massless Dirac fermions with vf = 3Jao/2h playing 
the role of the speed of light c, the quasimomenta K 
and are referred to as Dirac points [6, 32 . A lot of 
the excitement about graphene can be attributed to the 
promise of observing relativistic effects with solid state 
electrons. While for graphene c/300 the effective 

speed of light for atomic Dirac fermions in an optical 
lattice would typically be of the order of mm/s. 



IV. SINGLE PARTICLE DENSITY OF STATES 

We now turn to the fate of the Dirac points when a har- 
monic confining potential is added to the lattice. With 
the discrete translational symmetry broken, we can ex- 
pect to find both delocalized states with a well defined 
quasimomentum and localized states consisting of many 
quasimomentum components. It is therefore no longer 
meaningful to discuss the dispersion, and instead we look 
for evidence of the Dirac points in the single-particle den- 
sity of states (DOS) 



(5) 



Here the sum is over the eigenstates of the tight binding 
Hamiltonian 'HI '^/^n) = ^^nlV^n)- Numerically, we find p{E) 




FIG. 2: (Color online) The single-particle spectrum for an 
infinite honeycomb lattice. The two bands touch at the six 
Dirac points located at the edges of the first Brillouin zone, 
which is indicated by the hexagon in the plane E — —3 J. 



by binning the eigenvalues into small energy intervals of 
varying width. Counting the number of eigenstates in 
each interval gives a good approximation to the DOS in 
the middle of the intervals, provided the widths of the 
intervals are small enough to capture the variation of 
p{E) with energy, but large enough that fiuctuations are 
smeared out. 

Before investigating the DOS in the inhomogeneous 
lattice we first recall how the Dirac points are manifested 
in the form of the DOS in the absence of the trap. Since 
an analytic expression for p{E) exists for the infinite lat- 
tice, this also constitutes a test of the numerics. 



A. Homogeneous lattice 

For the homogeneous lattice the single-particle DOS 
per unit cell has the analytical form O [33] 



Po{E) = 



2 \E\ 



7r2 J2 ^ 




(6) 



where K{z) = /q^^^[1 — zsin^t] ^^'^dt is the complete 
elliptic integral of the first kind, and 



(i + ifi)' 



4|fl 
4|fl 



(1 + lfl)' 



\(E/J)^-lf 
4 



\{E/J)^-lf 



< 1 



1 < 



< 3 



lfl<l 
l<lfl<3 



(7) 



(8) 



The analytical DOS is plotted in Fig.|3] In the vicinity 
of the Dirac point {E = 0) the linear dispersion leads to 
a DOS which vanishes as po{E) oc \E\ with no band gap. 
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The van Hove singularities at £^ = ±J arise due to the 
saddle points in the single-particle spectrum at the edge 
of the Brillouin zone, halfway between neighboring Dirac 
points. We note that the DOS is symmetric around the 
Dirac point, po{—E) = po{E). The spectral symmetry is 
broken if next-nearest-neighbor tunneling is included in 
the tight binding Hamiltonian. 

The histogram in Fig. [S] is the numerically calculated 
DOS, which agrees with the analytical expression except 
for a large peak at £^ = 0. This additional peak is due 
to edge states, an artifact of our finite numerical grid. 
These zero energy modes are localized at the boundary 
of the system and appear because we confine the system 
in a cylindrical box. But they can be studied in graphene 
nanoribbons [6^ and could be constructed in an optical 
lattice by applying a repulsive potential at the edge of 
the cloud [T51. 




FIG. 3: (Color online) Single-particle density of states per 
unit cell for the homogeneous lattice. The solid line is the 
analytical expression for an infinite lattice, ([s]), while the his- 
togram is the binned density of states from the numerical 
calculation with rmax = 60ao, corresponding to 8792 lattice 
sites. The bin size is varied to resolve the details in the spec- 
trum. 



B. Inhomogeneous lattice 

We now turn to the combined lattice and harmonic 
trapping potential. In Fig. [4] we plot the binned density 
of states for a range of trap strengths. We make the 
following observations on the form of the DOS of the 
finite system: as the trap strength is increased from zero 
the characteristic valley around the Dirac point ai E = 
is gradually filled in, and the minimum is shifted to 
higher energies. For /^r^ax > 12J {tzal >3.3-10~^Jfor 
^max = 60ao) the valley has been replaced by a plateau, 
and as hc is increased further the length of this plateau is 
extended. The peak due to the edge states is shifted to 



E = ^^r^ax 5 expected for eigenstates localized at the 
edge of the cylindrical box. At the same time the peak is 
broadened due to mixing of the localized edge states with 
delocalized states in the same energy range. Lastly, the 
symmetry of the DOS is observed to be nearly conserved 
(apart from the edge state feature), but around an energy 
^0 = i^^^max > 0. such that p{Eo - E) = p{Eo + E). 
Below we explain each of these observations by analytic 
arguments. 



1. Low energy limit 

In the low energy limit, the lower band of the pure 
lattice has a dispersion resembling that of a free parti- 
cle with an effective mass, m* = fi^ {d'^ E^^/ dql\q^o)~^ = 
2h^ /ZJoq. Hence the low energy DOS is that of a 2D 
harmonic oscillator, p{E) = {E — £^min)/(^^*)^7 with a 
characteristic frequency cj* = ^/Hi/m* and a minimum 
energy £^min = —3J-\-hu* given by the infimum of the 
lattice spectrum offset by the zero point energy of the 
oscillator. The low energy DOS is therefore a linear func- 
tion of the energy for E > £^min- 



p{E) = 



~2J 



(9) 



2. High energy limit 

At high energies E ^ J the kinetic (and lattice) energy 
is negligible compared with the trap energy and Q re- 
duces to the potential energy of a 2D harmonic oscillator. 
The eigenstates of the trap potential energy operator are 
localized states with energies Ej = ^/^|rj| . These have 
been observed experimentally in a one-dimensional opti- 
cal lattice with harmonic confinement [26]. The DOS is 
then 



p{E) 



dN{E) dN{r) dr 



dE 



dr dE' 



(10) 



where N{E) is the number of quantum states with energy 
less than E and N{r) is the number of lattice points in 
a circle of radius r. Geometric considerations show that 
the DOS in the high energy limit approaches the constant 
value 



piE) = 



Sir 



3a/3 i^clq 



(11) 



Accordingly, ptva^ forms a plateau at Sn/SVS 4.84 
at high energies as affirmed by the numerical spectrum 
in Fig. [4] 

In the finite system the plateau in the DOS is observed 
to begin at £^ = 3 J and end at ^ = ^/^r^ax ~ Hence 
the plateau appears if /^r^ax > 12J. At higher energies 
the DOS decreases with increasing energy ultimately van- 
ishing at the largest eigenvalue in the spectrum, which 
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FIG. 4: (Color online) Binned density of states for four different trap strengths na^/ J = 5-10 ^,1-10 ^,1.5-10 ^ and 1-10 ^ 
(histograms). The solid line is calculated using a local density approximation for p{E) (see text). The dashed and the dotted 



lines are the low energy limit H and the high energy plateau (11 ), respectively, while the dash-dotted line represents the local 
density approximation for p{E)on a lattice confined in an infinite box. 



is approximately given by £^ma 
observations are explained in section [VT 



nicix ' ' T'hese 
below where we 



discuss an approximation to the spectrum based on the 
slow variation of the trapping potential on the scale of 
the lattice modulation. 



V. LOCAL DENSITY OF STATES 



ible as a large peak at r = rmax and E = \i^r'^Q^^- By 
our averaging procedure the van Hove singularities are 
rounded. This demonstrates that the Dirac physics of 
graphene is accessible in an inhomogeneous honeycomb 
lattice, provided local spectroscopic probes are available. 
For clarity the local density of states has been divided by 
47rrAr/(3\/3), which is the number of lattice sites in a 
radial shell between r and r + Ar. 



While the Dirac point in the global DOS is erased by 
adding a confining potential to the lattice we now in- 
vestigate if it survives locally by calculating the local 
density of states (LDOS), which is indicative of the local 
structure of the spectrum. Specifically, we calculate the 
angle- averaged LDOS 

p{E,r) = Y. r ^\i^n{rmE-E„), (12) 

by a binning procedure, where we add the probability 
densities of all eigenstates in a narrow interval of both E 
and r. 

As is clear from Fig.[5]the LDOS as a function of energy 
at a fixed r looks just like a local copy of the homogeneous 
lattice DOS displaced along the energy axis by the local 
harmonic potential energy ^/^r^ with the edge states vis- 



VI. LOCAL DENSITY APPROXIMATION 

If ttosc ^ Cio there is no appreciable change in the har- 
monic potential over several units cells, and an approx- 
imation where the lattice is taken to be locally homo- 
geneous can be expected to be good. With this and the 
suggestive form of p{E^ r) in mind we now construct a lo- 
cal density approximation (LDA) to gain further insight 
into the shape of the DOS. Semi-classically the local DOS 
for a unit cell at the distance r from the trap center is 
given by 

PLDA{E,r) = Y,5{E-E^{r)), (13) 
q 

where E^{r) = E^^ \nr'^. This is just the DOS for 
the homogeneous lattice shifted by the local harmonic 
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FIG. 5: (Color online) The local density of states as a function 
of energy and distance from the center of the trap for naQ = 
O.OIJ. For clarity the increase in p{E, r) due to the expanding 
number of lattice sites in the enclosed area as r increases has 
been removed. The large peak at r = rmax is due to the edge 
states. The harmonic oscillator potential energy is indicated 
by the dashed line. The analytical DOS for a homogeneous 
lattice is indicated at r = 0. 



potential energy, i.e. pLDA{E^r) = po{E — ^/^r^). In the 
LDA the global DOS is found by integrating Plda{E^ r) 
over the entire lattice, weighted by the number of lattice 
sites at each distance r. This is proportional to 27rr, and 
the DOS may therefore be approximated by 



Plda{E) 



^ Jui 



po {u) du. 



(14) 



The normalization constant M is chosen such that 
PLDA{E)dE gives the total number of lattice sites 
inside the radius rmax, and we have substituted u = 
E — ^Hir'^. The finite support of the homogeneous 
lattice DOS implies the lower and upper limits ui = 
max(— 3J, ^ — |/^r^ax) ^2 = min(3J, respec- 
tively. 

This LDA is plotted in Fig. [4] and shows a remark- 
able agreement with the numerically calculated DOS, 
apart from the edge states, which are not captured by 
the semi-classical estimate. The efficacy of the LDA was 
demonstrated for a three-dimensional cubic lattice with 
a harmonic confining potential in [34 , and the method 
should be valid in any optical lattice potential as long as 
the condition aosc ^ clq is satisfied. 

Based on the semi-classical estimate we can explain 
the following features of the DOS: 

Scaling: in (14) the integral only depends on the 



strength of the trapping potential through the lower limit 
ui. Hence if £^ < ^K.r'^g,^ — 3 J the value of the integral 
is only a function of E. This implies a universal form 
of i^Plda{E) in the limit of an infinite lattice such that 
pLDA for all energies. This agrees with ^ (in 

the limit where na^ <C J such that the LDA is valid) 
and with (11). The universal form of i^Plda{E) for an 



infinite lattice is indicated by the dashed-dotted line in 
Fig. [41 

Limits: Plda vanishes for E < —3 J, consistent with 
the analytical low energy estimate £^min above, when the 
zero-point energy of the trap can be neglected. The semi- 
classical DOS also vanishes at energies E > £^max- 

Plateau: the high energy plateau is also characterized 



by considering the limits in (14). If 3 J < £^ < i^r 



3 J the integral is over the entire homogeneous lattice 
DOS and equals a constant independent of E. Therefore 
Plda{E) = (271 / n) X const, in that case. This explains 
the beginning and the end of the plateau. The condition 
for the plateau to appear is 3 J < ^/^r^ax~3J or /^r^ax > 
12J. 

Symmetry: within the LDA we can understand the 
symmetry of the DOS as follows: if we neglect the small 
zero point energy hw* the center of the spectrum is given 
by Eq = (£^max + £^min)/2 = |/^r^ax- ^y another change 



of variable to v = Eq ± E - 
then be written as 



the DOS at ± ^ can 



Plda{Eo ±E) = / po{v)dv. (15) 

K, J-Eo±E 



By partitioning the integration interval the integral can 
be split into two part Plda = Plda + p¥da^ where the 
first part 

277- A/" pEo-E 

pinAiEo ±E) = / po{v)dv. (16) 

^ J-Eo+E 

is the same for both arguments. For simplicity we con- 
sider only £^ > 0. The second part is 

pi'nAiEo ±E) = / po{v)dv. (17) 

1^ J±Eo-E 

Since the homogeneous lattice DOS is symmetric about 
zero energy, po{—E) = po{E)^ it follows that Plda{Eo — 
E)=plda{Eo^E). 

It is important to stress that the symmetry of the DOS 
for the trapped system is a finite size effect. The same 
applies for the critical value of for the onset of the high 
energy plateau as well as the finite length of the plateau 
as a function of energy for a fixed trap strength. For 
an unbounded system the high energy plateau stretches 
to infinitely high energies and i^p{E) follows a universal 
form as discussed above. This is shown by the dashed- 
dotted line in Fig. |4j which represent Plda{E) in the 
limit where r^ax ^ y2(£^f3J)//^, such that finite size 
effects are irrelevant for the energies shown (note that 
E <C Vq /2hzaQ — 3 J is needed fo r the tight binding model 
to be applicable, c.f. Section III). It is worth noting 
that while the DOS for the finite system develops gradu- 
ally from that of the homogeneous lattice as the trap- 
ping strength is increased from zero, the DOS of the 
trapped, unbounded system is qualitatively different from 
its translationally invariant counterpart, owing to the di- 
vergence of the harmonic oscillator potential as r ^ oo. 
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This dramatic difference between the infinite system DOS 
for K: = and in the hmit k was also noted by Hooley 
and Quintanilla for a cubic lattice [23^. 



VII. FERMIONIC DENSITY PROFILE 

Above we have accounted for the spectrum of a single 
atom in honeycomb lattice with harmonic confinement. 
We have found that the Dirac points of the homogeneous 
graphene spectrum survive locally in the presence of the 
harmonic trapping potential. In this section we consider 
how this can be confirmed experimentally. While Bragg 
scattering has been applied with great success as a spec- 
troscopic probe of atomic quantum gases j35H37] , a cal- 
culation of the response of a many-body system to this 
kind of perturbation is beyond the scope of this work. 
Instead we look for evidence of the underlying relativis- 
tic physics in the density profile of a trapped gas, since 
this observable is universally available in experiments. 

For simplicity we concentrate on the density n(r) of 
a zero temperature, noninter acting Fermi gas with N 
atoms, since this only entails summing over the prob- 
ability distribution of the N lowest eigenstates 



N 



(18) 



An ideal Fermi gas is realized with a degenerate single- 
component (fully polarized) gas of ultracold fermionic 
atoms due to the suppression of p-wave collisions and the 
symmetry requirements imposed on the wavefunction of 
identical fermions by the Pauli principle. In Fig. [6] we 
plot the density at each lattice point as a function of 
the distance from the center of the trap. The density at 
distance r from the trap center can also be written as 



n(r) = f \{E,r)dE, 

J — OO 



(19) 



where the Fermi energy is fixed by the constraint N = 
f n{r)(Pr. In the center of the trap a band insulator with 
unit-filling is formed at sufficiently high particle number. 
By comparing with Fig.lSlone sees that unit-filling at site 



i requires > 3J-\-htvrf such that the integral in ( 19 ) is 



over the full DOS of the homogeneous lattice (displaced 
by the local oscillator energy). 

Based on a local density approximation for the 
fermionic density profile it has previously been suggested 
that the Dirac points emerge as a shoulder in the den- 
sity at a radius corresponding to half- filling [7 . This can 
be understood as the position in the trap, where the lo- 
cal Fermi energy £^f(^) = — ^A^r^ crosses the Dirac 
point located at zero energy in the homogeneous spec- 
trum, such that the integral in (19) covers exactly half of 



the displaced homogeneous lattice DOS. This prediction 
is confirmed by our calculation using the single-particle 
eigenstates of the tight binding Hamiltonian. The density 



profiles plotted in Fig. [6] show the anticipated shoulder 
at half-filling. 



0.75 




0.25 



FIG. 6: Density profiles of noninteracting fermions in the 
combined lattice and trapping potential at zero temperature. 
From the left to the right N = 100, 500, 1000, 2000, 3000, and 
4000. The trap strength is /cag = O.OIJ. 



VIII. CONCLUSION 

We have shown how a confining potential alters the 
spectrum of a single atom in a honeycomb lattice. Even 
though the eigenvalues of the tight binding Hamiltonian 
are significantly modified by increasing the strength of 
the trapping potential, the characteristic spectrum of the 
homogeneous honeycomb lattice survives locally in the 
trap, provided the confining potential varies over a length 
scale much larger than the extent of a unit cell. This 
means that it should be possible to observe graphene-like 
physics with cold atoms in a honeycomb optical lattice, 
and hence that this system can be used to implement a 
relativistic quantum simulator. 

We have studied the density profile of a single- 
component Fermi gas and shown that the Dirac points 
emerge as a shoulder at half-filling. In addition, the local 
density of states suggests that the massless Dirac quasi- 
particles can be directly manipulated by a local spec- 
troscopic probe. However, additional calculations are 
needed to conclusively demonstrate that the local dy- 
namics is governed by the Dirac equation. 

The single-particle density of states was fully described 
by a combination of analytical and semi-classical argu- 
ments. Importantly, the numerically calculated spectrum 
was reproduced with striking accuracy by a local density 
approximation based on the density of states of the ho- 
mogeneous honeycomb lattice. This implies that statisti- 
cal mechanics calculations of many-body systems in the 
combined trap and lattice potential can be done with- 
out resorting to numerical diagonalization of the tight 
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binding Hamiltonian, provided the trapping potential is ish Natural Science Research Council, 
slowly varying over the size of a unit cell. 
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